Quantum Monte Carlo

Levy constrained search in VMC?

2009-08-27T00:19:00.000-05:00

The Levy constrained search formalism is a useful way to reason about density functional theory. It finds the minimum energy using two nested searches. The inner search seeks the minimum energy wavefunction for a particular density. The outer search seeks the density with the minimum energy.

This could be directly implemented using VMC. The main ingredient is energy minimization subject to a density constraint. Has anyone done this?

It wouldn't be very efficient compared to optimizing the VMC energy directly, but it would make an interesting connection between the methods.

Why scientific computations fail

2008-12-14T14:49:00.002-06:00

(I was busy this fall with another project. Now it's time to return to scientific programming)

The previous post covered several levels of tests for assessing the correctness of scientific computation. Now let us look at reasons why the results might not match what is expected.

Programming errors
Mistakes in derivation and transcription of equations
Violations of assumptions or conventions
Errors from floating point
Errors and undesirable behavior from the algorithm (often problem-dependent)
- Algorithms that have been well-studied by numerical analysts
- Those that have not

Programming errors

The usual sorts of problems in programs (array indices, memory problems, etc). This category of problems can be tested with unit tests and can make use of tools and practices from mainstream software engineering.

Mistakes in derivation and transcription of equations

This is similar to the implementation not matching the specification in traditional software engineering. With scientific computation, though, the "specification" could be more precise in principle (equations and algorithms), and computer algebra systems could aid with performing and verifying steps in the derivation.

Violations of assumptions or conventions

Usually happens during derivation of the equations. A stronger connection between the derivation and the final code might help this.

Errors and instabilities from floating point

Accumulated round-off, loss of precision from subtracting nearly equal numbers, and other problems. For more subtleties than you every wanted to know about, visit W. Kahan's site.
One of the most annoying parts about floating point is lack of good model for understanding (and analyzing) floating point errors. I will return to this issue in a future post.

Errors and instabilities from the algorithm.

This category might be better divided into two subcategories - algorithms that have been well-studied by numerical analysts (like matrix operations, certain ODE's and PDE's, optimization), and those that have not. Many algorithms have been well-studied and have nice packages available to use as black-boxes. However, most scientific programs are a composition of well-known algorithms along with some problem-dependent features as well.

Often times the behavior of the algorithm depends on the details of the physical system (ie, is it a metal or insulator)?

This category also includes science-based and numerical approximations.

Now, how do we go from test case failure to deciding which of these categories the problem falls into (and then, of course, to locating the precise problem)?

Unit testing will help with the first category. But not so much with the rest. Unit testing is helpful when you want to make sure things haven't changed. The hard part is getting it right the first time.
The different levels of tests in the previous post should give some insight into which category the problem might occur.

First I want to look at various models of floating point error and see if we can gain some understanding of the intrinsic uncertainty in a calculation. This seems to be a necessary prerequisite for answering the question of 'how close is close enough?', in order to decide whether a test has actually failed or not.

Assessing the correctness of scientific computation

2008-07-13T21:54:00.000-05:00

In the previous post I promised some thoughts on assessing the correctness of scientific computation.

Computational science programs are typically assessed with a series of test cases, ranging from order-of-magnitude estimates ("those numbers don't look right") to comparing with experimental data. I've attempted to divide these test cases into categories, where each category shares similar expectations as to what 'close enough' means and has similar implications into potential causes for discrepancies.

Order of magnitude estimates from general principles, back of the envelope calculations, or experience with similar problems.
Analytically solvable cases.
Small cases that are accurately solvable
For small systems, or simple cases, there are often multiple solution methods, and some of them are very accurate (or exact). The methods can be compared for these small systems. An important feature of this category is any approximations can be controlled and the effects made arbitrarily small.
Results from the same (or similar) algorithm
Comparing the same system with the same algorithm should yield close results, but now there is additional uncertainty in the exact implementation, and the sensitivity to input precision (particularly when comparing results from a journal article)
Results from a different algorithm
This looks similar to the second category (comparing with exact or more precise algorithms), but this is the case where there a multiple methods (with different approximations) for handling larger systems. There are likely to be more approximations (and probably uncontrolled approximations) involved.

(In electronic structure there is Hartree-Fock(HF), several post-HF methods, Density Functional Theory, and Quantum Monte Carlo. And possibly a few others.) Now one has to deal with the all the above possibilities for two programs, rather than just one.
Experimental results
In some ways this is quite similar to the the previous category, except determining "sufficiently" close requires some understanding of the experimental uncertainties and corrections in addition to possible program errors.

The testing process for each case involves running the test case, deciding whether the program result is 'close enough'. If it is, proceed to the next test case. If it is not, work to discover the cause of the discrepancy. Next post I'll look at a list of possible causes.

Engineering scientific software

2008-07-12T01:20:00.002-05:00

Some pointers from Greg Wilson on engineering scientific software.

The first pointer is to the "First International Workshop on Software Engineering for Computational Science and Engineering".

A couple of Wilson's short summaries of papers from this workshop
(from recent research reading)

Diane Kelly and Rebecca Sanders: “Assessing the Quality of Scientific Software” (SE-CSE workshop, 2008). A short, useful summary of scientists’ quality assurance practices. As many people have noted, most commercial tools aren’t helpful when (a) cosmetic rearrangement of the code changes the output, and (b) the whole reason you’re writing a program is that you can’t figure out what the answer ought to be.

Judith Segal: “Models of Scientific Software Development”. Segal has spent the last few years doing field studies of scientists building software, and based on those has developed a model of their processes that is distinct in several ways from both classical and agile models.

What I found interesting:

Kelly and Sanders suggest that the split between verification and validation is not useful, and propose 'assessment' be used as an umbrella term. (I will return to the point in a later post)

The second pointer is to the presentations at TACC Scientific Software Days.

One of the presentations there was by Greg Wilson himself on "HPC Considered Harmful". There is one slide on testing floating point code (slide 15). In particular the bullet point testing against a tolerance is "superstition, not science". (And he's raised this point elsewhere)

This statement has bothered me since I first read it, and so the next few posts will explore some of context surrounding this question. The first post will list different methods for assessing the correctness of a program, the next post will categorize possible errors, and finally a post looking at 'close enough'.

Computing an unbiased inverse

2008-05-29T22:23:00.000-05:00

In a previous entry, Reweighting is biased, I had discussed the problem of a finite-sample size bias in reweighting. In this entry I will investigate a possible solution.

But first, I had asked if this was known in the literature. A comment by Lee Warren pointed to this article, Population size bias in descendant-weighted diffusion quantum Monte Carlo simulations, which discusses this very issue (in a slightly different context, but the root issue is the same). In that context the problem is even more severe since the straightforward solution (increase N) carries a heavier computational cost.

One of the crucial issues is finding a good estimator for the inverse (1/D). The naive estimator (simply do the division) is biased. So the question is: are there any other ways of performing division that could be used? Yes, in digital arithmetic there are a couple of algorithms for division using only multiplication and addition - Newton-Raphson and Goldschmidt (see the Wikipedia entry on digital division)

In general, to have an unbiased estimator, each sample of the random variable must enter the estimator linearly (alternately, we can say each sample should only be used once and then discarded)

The Newton-Raphson iteration for the inverse (1/D) is x_i+1 = x_i(2-D x_i). To make this unbiased, in each iteration the value of D needs to be an independent sample, and the values of x_i from previous iterations also need to be computed from independent samples. One can think of the algorithm forming a tree, with each iteration branching to two children (each child is one of the x_i computed by the previous iteration), until reaching the leaves representing the initial iteration (x₀ = T₁ + T₂ D)

I tested this algorithm, and it appears to give unbiased results. There is a catch, though. The result has a large variance. Using the average of all the samples with the naive estimator yields a bias that is much smaller than the variance of the unbiased estimator (hence the mean squared error is a more appropriate metric for the combination of variance and bias)

Open questions/issues:

Is there some modification of the algorithm that will produce a smaller MSE than the naive estimator? Or even an MSE that isn't too much larger?

How many iterations are necessary - how to analyze the convergence of this stochastic Newton-Raphson iteration?

Try out the Goldschmidt algorithm.

Reproducible Research

2008-04-06T00:56:00.000-05:00

There were a couple of recent posts on Greg Wilsons's blog about reproducible research. The first is a pointer to a paper about WaveLab. The paper is from 1995, although it appears the philosophy hasn't changed much looking at a more recent (2005) introduction on the website. The second post contains a reference to the Madagascar project.

These projects follow the work of Jon Claerbout on Reproducible Research.

Thoughts

This work is all about designing a workflow for computational research, and creating tools to support that workflow. Much of it seems quite banal, but part of the whole purpose of doing this is to free you from having to think about the boring and repetitious steps.
Levels of reproducibility
1. Yourself, right away (when tweaking figures, for instance)
2. Another research group member, or yourself 6 months later
3. Someone in another group, (a) able to run the code and get the same results, and (b) able to understand the code
4. Someone in another group, able to use your description to write their own code that gets the same results
The last level is the gold standard in scientific reproducibility. Experiments are reproduced from published descriptions (ideally, although the descriptions aren't always complete). Or theorists rederive intermediate steps leading to a paper's conclusion.

In software development, levels 1-3 are what counts. (It's not very useful software development if others can't run your code and get expected results.). The last level is similar to multiple implementations of a standard (eg, C++ compilers from multiple vendors.)

The referenced works seem to focus mostly on processing of raw data (images, seismic imaging, etc) and the processing steps to produce a final result.

How would this apply in other fields? Say, ones with smaller input data sets, but much higher computational demands.

The Madagascar project, for example, uses SCons to coordinate the steps in producing research output. What happens when one of the steps involves several weeks of supercomputer time? Two immediate problems that spring to mind - one probably doesn't want this step launched without warning by the script. Furthermore, the steps for launching such a job are likely to be very site-specific, making it difficult for others to reproduce elsewhere.

Somewhat related entries of mine: Keeping notes and Solving Science Problems

Notes from PyCon 2008

2008-03-19T22:50:00.001-05:00

I attended PyCon 2008 this past weekend, and had a good time meeting people in the Python commmunity.

Some notes

I tried to give a lightning (5 minute) talk about sympy and code generation. Unfortunately, Ubuntu 7.10 on my laptop (Dell Vostro) failed to display to the external output, so I was unable to present to the larger audience. However, I was able to give the talk to the people at the science BOF session. The slides are available here

At the science BOF session I met Michael Tobis, the author an interesting system called PyNSol. It uses python to describe the differential equations, and generates Fortran for performance.

The python 2 to 3 converter looks interesting. Not so much for the conversion of python versions, but for the infrastructure developed for the tool.

Enthought was showing off some multi-touch displays. Being able to interact with 3D visualizations (and Google Earth) via gestures for zoom, pan, and rotation made for a nice demo. The real question is if this type of display would be useful for actual use.

Code generation now uses templates

2008-03-17T20:34:00.000-05:00

In Quameon, the derivatives for the Gaussian basis functions are computed with sympy and the python code is generated. Creating the code that surrounds the generated code is a bit of a tedious procedure - the syntax tree needs to be constructed by hand. But no more! The surrounding code can now written in python and a parser (written using pyparsing) will create the syntax tree for you.

Some details:

Identifiers surrounded by percent signs (eg %name%) get turned into a template node type for subsequent replacement.

Example template:


class %orb_name%:
  def compute_value(self,alpha,v,r2):
    x=v[0]
    y=v[1]
    z=v[2]
    %value%
    return val

The code replaces %orb_name% with the orbital name (s,px,py,px,etc) and %value% gets replaced with the expression for the value of the orbital (the methods for the derivatives not shown - see codegen/primitive_gaussian/deriv.py in SVN for the full template.)

Why python?

2007-12-21T08:31:00.000-06:00

The author of the first comment on the Quameon entry wonders if Python will be too slow.

One reason to use Python is ease of programming. Hopefully the advantages of programming in Python will outweigh the performance penalty.

However, at some point, runtime speed will be an issue. Here's a list of steps to improve runtime performance (in increasing order of effort/expected return)

Use psyco. Put "import psyco; pysco.full()" at the beginning of your program and get an instant speed increase.
Use the Shed-Skin Python to C++ compiler. (I haven't tried this yet - don't know if it works with Quameon.)
Profile and rewrite slow parts in C/C++/Fortran. I've done this with the inter-particle distance computation in a classical MC code in Python, and gained some performance.
Code generation. Right now, the code for the Gaussian basis functions is auto-generated (no more computing derivatives by hand!!). I would like to move more of the basic formulas to this technique. Then it should be easier to retarget the code generation to another language (in theory, I haven't tried it yet).

Other tidbits:

Monte Python uses Python and C++ for QMC (paper, code). They wrote the time-consuming parts in C++, and use Python to glue it all together. The algorithm variants they tested (parallel distribution strategies for DMC) were achieved through changing the Python part of the code.

Code generation could potentially result in *faster* code than writing a general framework in C++ or Fortran. This is because the code can be specialized for a given physical system, exposing more optimization opportunities. For an example of code generation, look at the development SVN area of PyQuante - there is a program that will produce a C++ program to compute the energy for a molecule given in the input file.

Quameon

2007-12-05T23:32:00.001-06:00

I've been working on a Quantum Monte Carlo code implemented in Python. It's called "Quameon" (the name is a mixture of letters from "Quantum Monte Carlo in Python"). There's a SourceForge project for it - here are links to the SourceForge project page, the project web page, and the project wiki.

The code doesn't do much currently, but I mention the project now so I can write blog entries about various aspects of getting the code working.

I'm currently trying to validate that the code works by using HF single particle orbitals (obtained from PyQuante) and no jastrow factor - this should reproduce the HF energy. It seems to work for a few atoms (He, Be, C), after sorting out some normalization issues with the basis sets (apparently all basis sets are not normalized - this doesn't affect the energy, but does affect the eigenvector coefficients). I will need longer runs to be sure. The next step is to profile the code to look for any quick and easy performance tuning opportunities.

The first target (after getting the code running reliably) is to implement several forms for jastrow factors and several VMC parameter optimization algorithms.

Keeping notes

2007-10-04T22:38:00.000-05:00

Keeping a record of work is important for research. I use two main tools for this purpose. First, a text file with entries in chronological order (one file per month, with names like "Sept2007"). I write the date and then start taking notes for that day in the file.

Pro - easy to use, chronological view easily available
Con - math hard to write well, ideas spread over several days hard to view topically.

The other tool I use is a wiki. It's useful to keep data in a structured, hierarchical format that's easy to edit.

Pro - ideas arranged topically, can do better looking math
Con - must work topically from the start, no good chronological view (Yes, I know you can view changes in order, but it's not as easy to reconstruct what I was thinking on a particular day.)

I want a single tool that does both of these - that can view the notes both chronologically and by topic.
For work-flow, I want the ability to enter ideas and notes free-form, and perform categorization and annotation at a separate time.

Another additional feature that would be nice is recording other actions on the system (at least by reference).
One example is check-ins to source control.

Another possible example - What if gnuplot were to record all the commands issued (and kept a copy of all the files used?) Graphs could be re-created and played back at a later date easily.

I looked around a little bit to see if there are other applications that might meet some of these requirements (or at least be an improvement over text files and a wiki). One class of app is a note taking application. I looked at Tomboy. It does have a Latex plugin (though I could not get it to work). Latex side-rant: It's great for making nice looking mathematics, but you wind up with "equations under glass". It's pretty to look at, but you can't manipulate it further. For mathematical content, I want something with more precise semantics (Content MathML, for example)

Another class of tool is mind mapping software. I looked at Vym (View your mind) briefly, but I don't think I will like this type of software. Viewing a graph of links is nice, but not as the main working view - I prefer a wiki view (pages with links)

Black Box Reweighting

2007-08-19T21:56:00.000-05:00

In the last post, I wondered whether the black box reweighting method might help with the bias.

The Black Box Reweighting (BBRW) method computes weights using the observed distribution of samples, rather than the targeted distribution. (The author's website contains a link to supplementary material with further discussion.)

I ran a simple test, and the bias was slightly worse with the BBRW method, compared to standard reweighting.

The BBRW method does a fine job of solving the problem the authors intended (correct handling of a distribution that may be incompletely sampled, and lower variance of the result - the lower variance of the result was apparent in my tests), but it doesn't solve the reweighting bias issue.

The search for a method to correct reweighting bias continues...

Reweighting is biased

2007-08-14T22:52:00.000-05:00

The reweighting technique, where samples from one probability distribution are adjusted with a weight factor to compute averages from a different (but usually very similar) distribution, is biased for a finite number of samples.

The details can be found here.

Question for the audience: Is this known? I don't recall seeing this mentioned any time the method has been presented.

Some possible routes to fix the problem are listed in the paper, but alas, I haven't been able to turn any of them into a useful solution yet.

VMC optimization papers

2007-08-10T22:29:00.000-05:00

A recent comment suggested these papers on VMC optimization:

Alleviation of the Fermion-sign problem by optimization of many-body wave functions

Wave function optimization in the variational Monte Carlo method

I'd also add these to the pile of papers to read and understand:

Optimization of quantum Monte Carlo wave functions by energy minimization

A simple and efficient approach to the optimization of correlated wave functions

And since this post is alread about VMC optimization, some further questions:

What about simultaneous geometry and parameter optimization?

The focus seems to be on single geometries - what's likely more interesting in the future is a family of geometries - has anyone parameterized the VMC parameters (by the bond lengths or nuclear coordinates)?

Correct estimators in DMC without forward walking?

2007-08-10T21:09:00.000-05:00

Computing properties other than the energy using Diffusion Monte Carlo requires forward walking to get correct answers. That may change, if this paper - Hellman-Feynman operator sampling in Diffusion Monte Carlo calculations - is correct.

Although as I look at the paper more, the technique doesn't look any easier than the forward walking implementation by Casulleras and Boronat (the paper does mention this similarity). It appears to use similar data in a different combination. I would be interested to see how the noise, bias, and stability characteristics compare.

The Hellman-Feynman theorem seems to be getting lots of attention in the QMC world recently. It's also been used to construct Zero-Variance Zero-Bias estimators (most recently example here).

Gravitational N-body and the Art of Computational Science

2007-08-02T22:35:00.000-05:00

A rather ambitious project to describe a gravitional N-body code and how to write it in Ruby: The Art of Computational Science

As Greg Wilson points out, something like this would be nice to have in Python.

Efficient Global Optimization

2007-03-17T23:51:00.000-05:00

The optimization method mentioned last post is an improvement on Efficient Global Optimization, presented in the paper Efficient Global Optimization of ExpensiveBlack-Box Functions (also available from on one of the author's websites)

Another paper listed in the references, A Taxonomy of Global Optimization Methods Based on Response Surfaces(by D. R. Jones, one of the authors on the EGO paper) is a nice overview of various optimization methods. It presents each method and gives ways in which it can be fooled into choosing the wrong point as the optimum. A fix is given, which then leads to the next method in the list.

One of the core ideas is regression based on Gaussian processes (kriging). To learn more about Gaussian processes, see www.gaussianprocess.org. The two Jones papers listed above also show how the GP based methods fit in with other basis set methods.

Useful paper for VMC optimization?

2007-03-02T17:44:00.000-06:00

I ran across this paper, An informational approach to the global optimization of expensive-to-evaluate functions, and wondered if it might be useful for VMC optimization?

The type of problem being solved seems to fit:

Function is expensive to evaluate.
Function evaluation may be noisy.

Many optimization methods estimate an optimal point and evaluate the function at that point. The method described in the paper estimates the uncertainty in the knowledge of the function, and evaluates the function where it will best improve our knowledge of the function.

Some potential issues:

This method doesn't make use of gradient information, which may make it less competitive than methods that do use gradient information.

How well does it scale with the number of parameters? The authors present a method to keep the cost under control as the dimension increases, but is the method still effective then?

March meeting

2007-02-23T08:21:00.000-06:00

This year's APS March meeting has a couple of sessions on QMC: H4: Recent Advances in quantum Monte Carlo Simulations and V21: General Theory: Computational Quantum Monte Carlo Methods.

Another session of interest (for me, anyway) is A23: Focus Session: High Pressure I - Earth and Planetary Materials

I will not be attending the March meeting this year, as my wife and I are expecting our first child in mid-March.

Science Blogging Conference

2007-01-18T08:24:00.000-06:00

The North Caroline Science Blogging Conference starts this weekend. It looks like it covers some introductory how-to material, and then dives into issues such as promoting public understanding of science, teaching, and how does blogging interact with research.

If you're put off by the term 'blog' being overused and overhyped, think of the conference as discussing the potential and challenges of a (relatively) new one-to-many and many-to-many communication mechanism and how that could affect the interactions between scientists and the public, and the interactions between scientists.

Basis sets galore

2007-01-17T21:52:00.000-06:00

More basis sets that you can shake a stick at over at the EMSL Basis Set Exchange. (Not sure why you would want to shake a stick at basis sets, though.)

Schrodinger equation in the mass media

2006-10-26T08:08:00.000-05:00

The Schrodinger equation appears in the background in the second half of Weird Al's video White & Nerdy.

The video also shows rolling dice (briefly), if one wishes to stretch its relevance to include Monte Carlo as well.

The Other QMC

2006-10-20T21:12:00.000-05:00

Quasi-Monte Carlo is the other method commonly referred to by the QMC acronym (I will abbreviate it QuasiMC to minimize confusion). It's not really Monte Carlo since the sequence of points is not random, but the point sets do have a Monte Carlo-ish feel about them.

And in fact, it does share the property of arbitrary termination with Monte Carlo integration (consider a grid - you have have use all the points at a particular spacing. Stopping half way across the region of integration is not going to yield good results).

The reason people are interested in QuasiMC is that its convergence is better than the 1/sqrt(N) of Monte Carlo. It does this by using point sets that are more evenly distributed than random points - there are fewer clumps of points or large gaps that appear in random sequences (Search for 'Low-discrepancy sequences'. Here's the Wikipedia entries for a technical discussion and some sequences).

QuasiMC point sets are usually generated from sequences where the low-discrepancy condition can be verified theoretically, but it is possible to use the intuitive approach of making the points maximally spread out (see this presentation for an example).

The following algorithm will generate a set of N points that gives better convergence than random (at least it worked on a simple integrand in one dimension.)

compute some random points (>> N points)

pick a starting point and put it in the set S.

find the point that is the furthest from all existing points in S, and place that point in S.

repeat the step 3 until N points are selected.

(I'm not saying this an efficient way to generate the points, just that it can be done.)

For addition information, see also chapter 7.7 in Numerical Recipes

QMC wiki

2006-09-01T20:21:00.000-05:00

Check out the QMC wiki at http://www.qmcwiki.org

It looks like a promising resource for the QMC research community.

Optimal histogram bin width

2006-06-27T01:15:00.000-05:00

Kevin Knuth wrote a paper about finding the optimal number of bins to represent data in a histogram (Optimal Data-Based Binning for Histograms). He starts from a piecewise constant density model and finds the (Bayesian) posterior probability from this model (equation 36, which is actually the log of the posterior). The posterior function is then maximized to find the number of bins that best models the data.

The article also investigates the number of data points for a reliable estimation of the density. The recommendation is 100-150 points, if the distribution is Gaussian.

It would be interesting to apply this method to radial distribution functions. However the assumption of a constant volume for each bin is not met in this case. There are several ways this could be adjusted, but I'm not sure they are valid (scale each bin count by the volume, or use non-uniform bin spacing to maintain constant volume)

Alternately, the discussion references other algorithms for dealing with variable bin-width models (which may be better for resolving multiple peaks anyway).